On the role of confinement on solidification in pure materials and binary alloys 
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We use a phase-field model to study the effect of confinement on dendritic growth, in a pure 
material solidifying in an undercooled melt, and in the directional solidification of a dilute binary 
alloy. Specifically, we observe the effect of varying the vertical domain extent [5) on tip selection, by 
^ ' quantifying the dendrite tip velocity and curvature as a function of 5, and other process parameters. 

^ , As (5 decreases, we find that the operating state of the dendrite tips becomes significantly affected 

by the presence of finite boundaries. For particular boundary conditions, we observe a switching 
of the growth state from 3-D to 2-D at very small 5, in both the pure material and alloy. We 
demonstrate that results from the alloy model compare favorably with those from an experimental 
study investigating this effect. 
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I. INTRODUCTION 

. Dendrites are one of the basic microstructural patterns seen in solidified metals. The mechanical behavior of the 
solidified product is often decided by the length scales set by these patterns. Study of dendritic growth is therefore 
r , motivated by the need to predict these length scales. The fundamental quantities that completely describe the 
' ^ ' growth of a dendrite at steady state under a given set of external conditions are its tip velocity and radius, which 
together define the so called "operating state" . Despite considerable advances in the understanding of solidification 
science, discrepancies still arise when one attempts to compare theoretical predictions of dendrite operating states 
with experimental observations. We find that some of these discrepancies derive from differences between the ideal 
conditions assumed in theoretical treatments, and those experienced by materials under actual experimental situations. 

The first theoretical treatment of the "free" dendrite growth problem was presented by Ivantsov jj . He considered 
a pure dendrite, modeled as a paraboloid of revolution with tip curvature pup, growing into an infinite undercooled 
■ melt with temperature T — > Too far from the advancing tip. The dendrite was assumed to be isothermal at the melting 
, temperature T^, and to be growing along its axis at constant velocity Vup in a shape-preserving way. Ivantsov found a 
' solution to the thermal transport problem, in which the dimensionless undercooling A = [Tm — Too)/{Lf/cp) = I{Pe) 
\ where Pe = pupVtip/2D is the Peclet number, D is the thermal diffusivity, and X is the Ivantsov function. The 
temperature has been scaled by the characteristic temperature Lf/cp, with Lf being the latent heat of fusion and Cp 
the specific heat. 

This solution presented a conundrum, because it showed that the transport problem alone did not uniquely specify 
the operating state of the dendrite, i.e., the single combination of pup and Vup observed in experiments. Additional 
I considerations, such as the effect of curvature on the melting point stability and eventually the anisotropy of 
'"^ the surface tension led to a second condition a* — 2doD / p^^pVup where a* is called the selection constant, and 

! rfo is the capillary length. The combination of Ivantsov's solution (modified for surface tension and its anisotropy) and 
O 1 the condition a* is constant gives a unique operating state. Numerical simulations using the phase-field method |8|| at 
J-^ ' large values of A, have found agreement with the predictions of this body of work, known as microscopic solvability 
^ theory. 

Glicksman and co-workers developed experimental techniques for studying the solidification of pure materials, with 
I^S the objective of observing the operating state. They performed experiments with phosphorous and transparent 
' analog alloys like succinonitrile (SCN) iu\ and pivalic acid (PVA) [ij. The results of these careful experiments found 
some areas of agreement with microscopic solvability theory, in particular, the value of a* was found to be constant, 
but the operating combination of pup and Vup did not agree. Provatas, et al. were able to explain this discrepancy by 
showing that for the low undercooling conditions found in the experiments, interaction between neighboring dendrite 
branches ^2>^3 affected the operating state. 

Experiments have also been performed to examine the role of superimposed fluid flow on dendritic growth. Gill 
and coworkers 0, used SCN in a special cylindrical chamber with a bellows to effect fluid flow. Bouissou and 
Pelce 13 performed experiments with a flowing alloy of PVA and a seed conflned between microscope slides. Saville 
and Beaghton |17| pres ented a theoretical analysis which extended Ivantsov's solution to consider the superimposed 
flow. Jeong et al. [l8| performed phase-fleld simulations of these experiments, and once again found discrepancies 
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with theory. They conjectured that the differences arose because of the effect of finite containers in the experiments, 
leading to boundary conditions which differed from the assumptions of infinite media used in the theory. 

Dendrite tip theories for constrained growth, such as directional solidification of dilute binary alloys between 
microscope slides, have been developed by Trivedi Il9j and, Kurz and Fisher 20j. They have shown that a relationship 
of the form plipVup = constant, should hold for constrained growth just as in free dendrite growth. Early experiments 
by Somboonsuk et al. |2lj | in samples with slide separation greater than 150 /im have shown excellent agreement with 
this theory. However, in recent studies Liu et al. l22f| have demonstrated that experimental results start to deviate 
significantly from theory when the slide separation approaches the scale of the primary dendrite spacing. 

In this article, we systematically study the role of confinement on dendritic growth. Since every experiment is 
performed in a finite container, we feel that this effect cannot be ignored. The first two authors have previously 
reported a study on confined growth in pure materials 23] . Here we extend our investigations to directionally 
solidified binary alloys. For purposes of continuity and completeness, we will again present results from our study on 
pure materials. For the binary alloy, recent experimental data from Liu et aL|22l|. will provide us with an avenue for 
testing our numerical predictions. 

II. MODELING 
A. Phase field model, adaptive grids and numerical methods 

The objective in a general solidification problem is to solve the equations governing thermal and solute transport, 
subject to boundary conditions on the solid-liquid interface (moving boundary) and elsewhere. If melt convection is 
to be modeled, one needs to solve the momentum equations for fluid flow simultaneously with the above transport 
equations. Imposing the interface boundary conditions upon discrctizing the governing equations poses a difhculty 
however, since the interface, as it evolves, will not in general align itself with a fixed set of mesh points. 

The phase-field method eliminates the sharp liquid-solid boundary by introducing evolution equations for a continu- 
ous order parameter cf) g [—1, 1], where (p = —1, -|-1, corresponds to liquid, solid and interface respectively. Thus, the 
arduous task of solving the transport equations separately in liquid and solid domains while simultaneously satisfying 
boundary conditions on arbitrarily shaped interfaces, is replaced by that of solving a system of coupled differential 
equations; one for the evolution of cf) and one for each of the transport variables (temperature, concentration and ve- 
locity). Phase-field modelin g ha s been an active area of research in the past decade, and we refer the interested reader 
to original work by Langer '24'| , Karma and Rappel , and Beckermann et al. [23 for derivations of the phase-field 
equations and selection of phase-field parameters ensuring convergence to the original sharp interface problem. 

The phase-field model introduces a parameter Wq that connotes the finite width of the now 'smeared' interface. 
Karma and Rappel showed that the model converges to the sharp interface equations when p = WoVup/D <^ 1, 
where D is the thermal or concentration diffusion coefficient, and Vup is the nominal tip velocity of the dendrite. 
Resolving the interface on a discrete mesh requires that the mesh spacing Aa; ~ Wq, while demanding that the 
diffusion field not interact with the boundaries leads to the domain size Lb 3> D/Vup- Satisfying these requirements 
causes calculations on regular meshes to quickly reach the limit of available computing resources. For example, if we 
choose p — 0.01, fix LbVup/D = 10, and enforce Ax — Wq, then we find that the number of grid points per dimension 
on a regular mesh should be at least Lb /Wo — 1000. This makes computations challenging on regular meshes even 
in 2-D, while 3-D computations may not be practical at all, depending on available computing power. 

We have mitigated this problem successfully by solving the equations on an adaptive finite element mesh 0, . 
In three dimensions, we use eight-noded trilinear brick elements stored using an octree data structure. A local 
error estimator indicates refinement or coarsening of the mesh, and this permits tracking of the interface as well as 
resolution of gradients in the other fields. There are six degrees of freedom at each node (three velocities, pressure, 
temperature/concentration and 0), and a typical computation reaches well over one million unknowns. The finest 
elements (Axmm), which are distributed near the interface, now need to be order of Wq. 

For our studies on pure materials, we have used a finite element discretization of the 3-D phase field model developed 
by Karma and Rappel 01 ■ In order to account for the effects of melt convection we adopt the formulation presented by 
Beckermann et al. |25| . who use an averaging method for the flow equations coupled to the phase-field. By appropriate 
choice of phase-field parameters we have ensured zero interface kinetics, which is a valid assumption for the range of 
undercooling we are concerned with. 

For our alloy simulations, we have used a one sided (vanishing solid diffusivity) phase-field model [26L l2^ ^28*1 . with 
a frozen temperature approximation. In a directional solidification arrangement, for certain values of the problem pa- 
rameters (particularly when simulating real materials), a considerable amount of time can elapse before the transients 
vanish and the solid-liquid interface reaches steady state. In that time, the interface can encounter the end of the 
simulation box if the equations are solved in a reference frame that is fixed globally, and if the box is not large enough 
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to contain the diffusion field. To alleviate this difficulty, we have solved the phase-field equations in a coordinate 
frame translating with the pulling speed. This saves some computational expense by allowing us the use of smaller 
boxes. We have not investigated the effect of melt convection in our numerical experiments with alloys. 

In a recent article, Echebarria et al. have emphasized that the same choice of phase-field parameters that 
produced zero interface kinetics in the pure material cannot also ensure this condition in the alloy model. This 
is due to the presence of certain additional terms in the kinetic parameter P, that arise out of accounting for the 
discontinuity in the concentration field at the interface in a model with vanishing solid diffusivity. To ensure that 
the kinetic coefficient is negligible at the interface, the phase-field relaxation time r needs to be made temperature 
dependent in this region by setting 



Here tq is the usual relaxation time, k is the partition coefficient, z is the distance from the interface, Vp is the pulling 
velocity, t is time and It is the thermal length p^ . 

We evolve the nonlinear order parameter equation using a Forward-Eulcr time stepping scheme, while the linear 
thermal/solute transport equations are solved using the Crank-Nicholson scheme with a diagonally preconditioned 
conjugate gradient solver. The transport equations typically converge in fewer than five iterations per time step. The 
3-D flow equations for the pure material are solved using the semi-implicit approximate projection method (SIAPM) 
|29j| . Details of the above numerical methods and the finite element formulation are omitted here, as they have been 
presented elsewhere [Tsj . 



Our three dimensional simulation domain is the rectangular parallelepiped illustrated in Fig. ^ with edge lengths 
along the x, y and z axes; L^;, Ly and S respectively. The edges of the box are oriented along (100) cubic crystallo- 
graphic directions. 



r = To l-(l-fc) 



z - V^t 



(1) 



B. Geometry, initial and boundary conditions 




FIG. 1: Simulation Domain. All surfaces are modeled as symmetry planes in pure material simulations, whereas the surfaces 
2/ = and y = Ly are periodic boundaries in the alloy simulations. 
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Table 1. Pure material simulations 



Case 


A 




S 


l(a)-(f) 


0.55 





128,64,32,16,8,4 


2(a)-(f) 


0.55 


5 


128,64,32,16,8,4 


3(a)-(e) 


0.25 





128,64,32,16,8 


4(a)-(f) 


0.25 


5 


128,64,32,16,8,4 


5(a)-(d) 


0.15 





128,64,32,16 


6(a)-(e) 


0.15 


5 


128,64,32,16,8 



For the pure material, the initial condition is a spherical solid seed with a radius greater than the critical nucleation 
radius, centered at the origin depicted in Fig. ^ Because of the inherent symmetry in the growth of the seed, it is 
usually sufficient to model one octant in three dimensional space. However, if forced fluid convection is incorporated 
along a particular direction, then the solidification rates into and counter to this direction become unequal. For 
example, if there is a flow parallel to the x axis, x = is no longer a plane of symmetry. To account for this break in 
symmetry, we have to model at least a quadrant of space, with z — and y = as planes of symmetry. 

In the simulations with the pure material, the dimensionless thermal field u is subjected to zero flux (Vu • n = 0) 
boundary conditions on all surfaces. Fluid flow, when it is included in our study, is imposed as an inlet boundary 
condition Uoo, normal to the face x — —Lx/2. The velocity field is subjected to symmetry boundary conditions on the 
domain walls, and is forced to vanish in the solid ((/) — 1) by an appropriate formulation of the momentum equations 
(see ll^)- We fix the lateral dimensions of the simulation box (L^ ~ 512 and Ly = 256 are typical values), and study 
the interface evolution as a function of (5, which is varied from 128 to 4. Here, ij,, Ly and 5 are in units of the interface 
width Wo- For very small S (< 8), steady growth conditions are reached relatively quickly for large undercooling. 
To save on computational cost in these runs (where Axmin = 0.5), we sometimes use shorter lengths for and Ly, 
chosen to ensure that the diffusion field does not interact with the ends (in the x and y directions) of the box. For 
smaller A however, it typically takes much longer to reach steady conditions, and when melt convection is included, 
it can take impractically long CPU times to get converged results. For these cases, we terminate our runs when the 
tip radius/velocity versus time curves start to even out. Fortunately, it turns out that the behavior we are interested 
in appears for combinations of S and A where steady state conditions are always achieved. 

In the alloy simulations, the initial condition is a planar interface at a: = Xq, perturbed by randomly spaced finite 
amplitude fluctuations. The box in these simulations is taken to represent the shallow channel between microscope 
slides where directional solidification conditions are imposed, viz. a fixed thermal gradient moving at a constant speed 
Vp. Once again, in this arrangement we study the influence of the depth of the channel S (or equivalently the sample 
film thickness), on interface morphology. To minimize the diffusion fleld's interaction with the lateral boundaries, Lx 
and Ly are chosen to be relatively large 256). We enforce zero flux boundary conditions on the concentration field, 
on the surfaces x — Lx/2 = —L^jl and z = = 5, while periodic boundary conditions are imposed on the boundaries 
y — ^ and y — Ly. The rationale behind periodic boundary conditions is to be able to simulate an infinite domain in 
V- 

Unless otherwise stated, on each boundary, we employ the same type of boundary condition on the phase-field 
variable </>, as we do on the transport variable. Where V0 • n = 0, the material "wets" the boundaries, and the 
corresponding contact angle is 90°. Semoroz et al. have previously used this technique to capture wetting of solid 
surfaces, with a two-dimensional phase- field model for binary alloys [s^l- We also show a calculation with (f> — —1 
on the boundaries, which is equivalent to making the material "non- wetting" (contact angle — 0°). The real contact 
condition probably lies somewhere in between these two extremes. 

III. EFFECT OF SMALL 5 IN A PURE MATERIAL 

In this section, we report the effect of changing S on the tip of a pure material dendrite, evolving along the negative 
X axis (upstream direction when fiow is present). We simulated the cases shown in Table 1. We used a fixed value 
for the four- fold anisotropy in all our simulations (54 = 0.05). We have not corrected for grid anisotropy Q in these 
calculations, but work at a grid spacing where its effect is known to be small [l^ . 

We use the following values for parameters in our calculations: interface width Wq = 1, time scale for interface kinetics 
To = 1, coupling constant A = 6.383, thermal diffusivity D ~ A, capillary length do = 0.13851Vo (which leads to zero 
interface kinetics), and Prandtl number Pr = 23.1, where Wq, tq and D are in dimensionless units (see 
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A. "Wetting" boundary conditions, • n = 



To make ideas more concrete, we choose cases 3 and 4 as a representative subset of our computations and present 
detailed analyses on those runs. Figures [21 and |21 show the upstream dendrite's tip velocity and radius respectively, as 
functions of the box height S. 
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FIG. 2: Tip velocity vs. S, corresponding to cases 3 and 4 in Table 1. 



In these plots, S and pup were made dimensionless by scaling them with do, and Vup is scaled by D/do. We compute 
the tip radii p^z and p^y along two principal planes using the method of Jeong et al. , and estimate the mean tip 
radius by the formula 

Pup = 2 (— + — ) . (2) 

\Pxy Pxz J 

We can see from Fig. |21that for large values of (5, the tip velocity remains relatively unaffected by the box height. 
However, as we go to very small heights Vup decreases quite dramatically. As 5 is decreased, there is also a gradual 
decrease in the tip radius pup- Clearly enough, box height has a pronounced effect on tip dynamics. Fluid flow 
induces a parallel shift in these curves. The dendrite tip velocity increases uniformly in the presence of flow |l8l | . On 
the other hand, the tip radius is lower than the case with pure diffusion. 

The observed trends can be explained as follows. As long as 5 is sufficiently large, the thermal field enveloping 
the dendrite will interact with the upper boundary at a distance that is relatively far behind the tip. In particular, 
the thickness of the thermal boundary layer near the tip remains unaffected by this interaction. However, as 5 is 
decreased, this thickness can grow quite rapidly. We illustrate this effect by examining the temperature profile in the 
x-z plane, as shown in Fig. ^ It is evident that the temperature contours are more spread out in Fig. |4(b)| where 
i5 = 8, compared to those in |4(a)| where 5 = 64. The increased boundary layer thickness, decreases the thermal 
gradient into the liquid at the liquid-solid interface, which in turn retards the growth rate as a direct consequence 
of the Stefan condition. Due to the zero flux boundary condition on the plane z = (5, further reduction in 5 makes 
heat transfer in the vertical direction almost completely ineffective. Tip curvature in the x — z plane vanishes and the 
dendrite switches morphology from 3-D to 2-D. We note that once the dendrite goes 2-D, pup = Pxy 

An interesting result here is the 3-D to 2-D transition. We have performed tests with finer meshes (more elements in 
the vertical direction) to ensure that it is not simply an artifact of poor grid resolution. We believe this phenomenon 
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FIG. 3: Tip radius vs. 5, corresponding to cases 3 and 4 in Table 1. 




(a) (S = 64 corresponding to Case 4(b). u ranges from A = -0.24 to L = 0.008 




(b) 5 = 8, corresponding to Case 4(e). u ranges from A = -0.24 to L = 0.008 



FIG. 4: Temperature contours for two different box heights when A — 0.25 and Uoo = 5.0. In each case, the letter X symbolizes 
the dendrite outline (bold contour), under steady growth conditions. Contours are plotted in intervals of 0.025. Notice that 
the contours near the dendrite tip are more spread out when 5 is smaller. 



7 




(a) <5 = 16, 3-D dendrite 



(b) 5 = 4, 2-D dendrite 



FIG. 5: 3-D to 2-D dendritic transition at small S. The shaded surface is the dendrite (0 = 0). For these runs A = 0.55 and 
(7oo =0. 



to be a consequence of the V(f> • n = boundary condition on the upper boundary. The only way for the solid-liquid 
interface to match this boundary condition at small S is for the curvature in the x-z plane to vanish. This is illustrated 
in Fig. 

It is also interesting to observe the effect that melt flow has on the interaction between the tip and the boundary. 
We find that melt convection reduces the value of S where the tip velocity deviates from its nominal value. A 
straightforward explanation for this effect is that advection increases the rate of heat transport from the upstream 
dendrite arm. This increases the growth rate (hence Vup), while compressing the boundary layer, whose thickness 
scales as D/Vup. Thus the dendrite remains three dimensional for a smaller value of S than was previously possible, 
with convection absent. A more negative value of the undercooling also has qualitatively the same effect on the 
strength of the tip-boundary interaction, since Vup again increases in this case. Fig. |S1 summarizes these observations 
succinctly. Both melt convection and larger undercooling, cause points on respective curves whereupon interaction 
effects become important, to shift to the left. 

Using a graph such as Fig. (BJ one can derive a semi-quantitative estimate as to when the operating state becomes 
affected by the finite height of the container. If one assumes that the tip velocity at (5 = 128 for each case is 
approximately the tip velocity of a dendrite growing under identical conditions in an infinite domain, it is possible 
to quantify the influence on the operating state in terms of a percentage deviation in the true tip velocity from this 
nominal value. If, for example, we consider deviations of the order of 3% to constitute a change in the operating state, 
a least squares fit to these cut-off points on the respective curves in Fig. yields the criterion 



If this condition is not satisfied, then it is likely that the operating state is being influenced by the boundary. A simple 
condition such as the one in Eqn. Q may be used as a rule of thumb in determining if experimental studies on free 
dendrite growth, in geometries similar to ours, are free from contamination. In fact, one may have intuitively guessed 
a condition of the type S > a (D/Vup) (where a is some constant) to apply based on physical arguments alone, and 
Eqn. (O supports this conjecture. 



To underscore the importance of phase field boundary conditions in selecting a particular growth state, we present 
results from another simulation with A = 0.55, Uoc — and 6 — 4. This time we impose (f> — —1 on the upper 
boundary, which corresponds to a physical situation where the solidifying material is not allowed to wet the surface. 
A three dimensional surface plot of a steadily growing dendrite is shown for this case in Fig. 13 

As expected, the dendrite does not adhere to the top surface at all. We conclude that this will be the case for 
any S if Dirichlet conditions of this nature are imposed. It is evident then, that the 3D-2D transition that we saw 
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B. Non-wetting boundary conditions 
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FIG. 6: Tip velocity as a function of & for different undercooling and flow conditions. A weak power law relationship emerges 
between the 5 below which interaction effects are strong, and diffusion length D/Vup. 




FIG. 7: Three dimensional dendritic growth for A = 0.55, Uoo = and 5 — A, with ( 
the dendrite does not wet this surface. 



-1 on the upper boundary. Note that 



previously is a strong function of the boundary conditions imposed on the phase field. It would be interesting to 
conduct a detailed investigation of the effect of different boundary conditions on the tip in 3-D. The interested reader 
is referred to the article by Semoroz et al. [s^l, which discusses the influence of different contact angles on a dendrite's 
tip velocity, in two-dimensional thin film solidification. 



IV. EFFECT OF SMALL 5 IN A DIRECTIONALLY SOLIDIFIED ALLOY 

In this section, we describe results from our simulations on a directionally solidified alloy, and make comparisons 
with the experimental data of Liu et al. ■ These simulations were conducted for the range of pulling speeds and 
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Table 2. Alloy simulations 



Case 




S 


1 


0.8 


64,32,16,8,4 


2 


1.0 


64,32,16,8,4 


3 


1.5 


64,32,16,8,4 


4 


2.0 


64,32,16,8,4 



Table 3. Physical properties of a SCN-Salol alloy system 



\rn\ (Liquidus slope) 0.7 K/wt. % 

D (Diffusion coefficient) 8 x 10~^° m^/s 



r (Gibbs-Thomson coefficient) 0.64 x 10~^ K m 

k (Partition coefficient) 0.2 

G (Thermal gradient) 4 K/mm 

do (Capillary length) 3.265 x 10"* m 

It (Thermal length) 4.9 x 10"* m 



channel depths shown in Table 2. The depth S was varied from 64 to 4, and in each case, the simulation was continued 
until the interface became stationary in the moving frame. From the stable array of cells formed, a few were chosen 
as representative of the array, their tip radii were extracted in the two principal planes using a least squares quadratic 
polynomial fit, and the mean tip radius of each cell was computed using Eqn. (j^J. An average over these values was 
taken to be the mean tip radius of the interface. We note that this process required some approximation, especially 
in cases where we found steadily growing cells of disparate sizes that would have made such averaging inappropriate. 
In such cases, our usual approach was to choose cells farthest from the boundaries, and where even this failed to 
provide clear-cut choices, smaller cells were chosen because of their lesser likelihood to split. In a majority of our runs 
however, the choices were unambiguous. For certain values of d we found interface evolution to occur in a way that 
cells would creep on either the top or bottom surfaces. In those cases, the phase-field boundary conditions on the 
respective surfaces allowed us to treat them as symmetry planes for the purpose of calculating tip curvature. 



A. Selection of simulation parameters 

The following values for the lateral dimensions were seen to yield satisfactory results. = Ly = 256 when S < 16, 
and Lx — 256, Ly — 128, otherwise. We chose our simulation parameters to keep computations tractable. It took 
about 90 hours of CPU time on a 3.1 GHz processor to simulate a typical directional solidification experiment for a 
chosen set of phase- field parameters on a mesh with about 170 000 elements {6 — 64). The interface required about 
250 dimensionless time units to reach steady state in this case. As noted earlier, the use of a moving reference frame 
allowed us to cut substantial costs associated with the need for larger domains to prevent the diffusion field from 
running out of the domain. 

We did not attempt to model a real material in this study as this caused our simulations to become considerably more 
expensive. To illustrate this, consider a SCN-Salol system having the properties listed in Table 3. The conditions in a 
directional solidification experiment are completely described by the following two dimensionless control parameters, 
M = do/lx = 6.66 X 10~^ and S — Vpdo/D = 2.04 x 10"'*, where do is the capillary length, It is the thermal length 
and D is the solute diffusivity in the liquid phase; for a pulling velocity of = 5/im/s. To get converged results 
with the phase- field model we require that the solution become independent of the parameter e — Wo/do- After 
ensuring vanishing interface kinetics, the following relationships involving the phase-field parameters are realized: 
Dto/Wq = aia2e, VpTo/Wo = Sa^a2e\ a,nd It/Wq = l/{eM). Here, oi = 0.8839 and 02 = 0.6267, are constants that 
arise in the phase-field formulation. [23 

Echebarria et al. have shown that mesh converged results can be obtained with e as large as 50. Setting 

e = 50, gives us an under-determined system of three equations with the five unknowns D, Vp, It, tq and Wq. Making 
arbitrary choices for two of these parameters by setting Wq = tq = 1, we obtain D — 27.7, Vp = 0.2825, and It — 300. 
A large value of It implies needs to be very large at steady state, even in a moving reference frame, to contain 
the diffusion field. To avoid this, if we choose a more tractable value for It (say 100), and fix tq = 1, we now get 
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Wo = 0.3333, D = 3.077 and Vp = 0.094. Thus, the smallest element in our mesh needs to be at the very least 
Ax = 0.3333. Stability considerations now place a severe restriction on the size of the time step (At) needed for 
solving the phase-field equations by the Forward-Euler method, as At ~ Ax'^ . Since it is clearly impossible to choose 
both It and Wq independently, calculations involving real materials are typically more expensive. 

We choose instead a more computationally favorable set of dimensionless parameters M and S for our study. To 
achieve our primary objective, which is to study the effect of S on interface morphology, we anticipate, and in the 
following paragraphs demonstrate, that this hypothetical treatment will not obscure any physics. The parameters 
used in our study are: tq = Wq = 1, D = 20, k = 0.8, £4 = 0.05, and It = \m\{l — k)Coo/kG = 50; where m is the 
liquidus slope, G is the imposed thermal gradient, and Coo is the far field solute concentration. The condition for 
negligible interface kinetics gives do = 0.0277 and therefore e = 36.1, which is sufficiently small to ensure convergence. 
The size of the smallest element in our adaptive mesh is Ax — 1, when S > 8, and Ax = 0.5, when S — A, while 
At = 0.005 is the size of the time-step. For these parameter choices, the dimensionless control parameters work out 
to be M = 5.54 x 10^^ and S - 1.385 x lO'^. 



B. Interface morphology and comparison with experiments 



In order to test the model, we initially performed a set of runs to verify the MuUins and Sekerka stability limit of 
a planar interface puj l. perturbed by small sinusoidal perturbations. We found that the model captures the stability 
spectrum correctly. Having convinced ourselves that this fundamental requirement was met, we proceeded with our 
study. Comparisons with the experimental data in this section offer a better validation of the model. 

Fig. IHl shows the computed interface morphology at different values of S. For small values of pulling speed (Vp < 2), 
the steady state consists of a stationary array of cells as in Fig. |S1 However, unlike the cells observed in experiments, 
that are usually characterized by blunt tips, these appear to have sharper and better defined tips, giving the impression 
of dendrites. It is conceivable that ignoring thermal noise in our calculations is responsible for the absence of side- 
branches on these structures, that are typical of dendrites. At large values of S, we find that the tip radii of these 
cells, measured on the two principal planes, are almost identical. However, as S decreases, the in plane radii diverge 
from one another. In particular, the radius in the x-z plane becomes significantly smaller, and cross sections of the 



cells look elliptic. At (5 = 4, for small pulling speeds (Vp < 1), we get a two-dimensional interface (Fig. 8(c) I. The 
inter-cellular spacing also increases as 5 is decreased. 

As pulling velocity is increased, the morphology becomes finer, with sharper and more tightly packed cells. This 
behavior is consistent with that seen of both cellular and dendritic arrays in directional solidification experiments 
[21, 22, 32], where the primary spacing decreases with Vp. 

For the set of phase-field parameters we have chosen, if > 2.5, the interface does not reach steady state in a 
reasonable amount of CPU time, due to repeated tip-splitting of the cells. Splitting is initiated by oscillations that 
appear at the tip and propagate downward along the trunk of the cell. Cell spacing and shapes change very rapidly 
in this regime. We did not continue these runs any further to check if steady state is reached eventually. Instead, 
we set Vp = 2.0 as an upper bound on the pulling velocity, below which a stationary state was always the outcome. 
Fortunately, this still left us with sufficient sample space to conduct our study and make effective comparisons. 

To enable plotting of our results on the same graph with the experimental data of Liu et ai, which corresponds to 
a SCN - 0.7 % wt. Salol system (properties in Table 3), we non-dimensionalized the axes as follows. The abscissa 
is the pulling speed (i.e. the tip velocity at steady state) Vup, scaled by a characteristic velocity D/dok, while the 
ordinate is the tip radius pup, scaled by the diffusion length D/Vup. One may appreciate the fact that the abscissa is 
in fact the dimensionless parameter that we had earlier denoted by S, multiplied by k. Fig. shows a comparison of 
the data. The open symbols correspond to the experimental data of Liu et ai, while the solid symbols correspond to 
our calculations. A comparison between our data and theirs holds up surprisingly well. Of special significance are the 
following two observations: 1) Although we conducted our simulations at values of S and M that were each about 
an order of magnitude off theirs, the two sets of data correlate very well, i.e. appear to collapse on parallel curves 
that are not significantly different by way of intercept. This tells us that our choices of parameters for scaling the 
axes are appropriate. 2) Since their experimental data correspond to dendritic arrays, the cell-like structures we have 
computed are likely branchless dendrites. 

In their experiments, Liu et al. note that tip radius data for dendritic arrays agree quite nicely with a relationship 
of the form PtipVup = C, where C is a constant dependent on do, k and D, as postulated by theoretical models of 
constrained growth However, when S is of the order of inter-dendritic spacing Ai, this agreement deteriorates. 

This is evident in Fig. |3 where the p^^pVtip = C line is shown as a guide to the eye. The open circles, which are 
data for S = 12.5 /im deviate in both slope and intercept from this line, which passes through the rest of their data, 
indicating a breakdown in the relationship. We observe similar trends in our data, viz., the line pl^pVup — C fits our 
data at (5 = 32 and 64 reasonably well, but as 6 decreases from 16 to 4, this agreement deteriorates. 



11 




(b) <5 = 16, Vp = 1.5, three-dimensional cells 




(c) (5 = 4, Vp = 0.8, two-dimensional cells 

FIG. 8: Interface morphology in a directionally solidified alloy. At steady state, stable arrays of three-dimensional cells appear. 
Note that at 5 = 4, the array comprises of two-dimensional cells. 



As in the numerical experiments with the pure material, we find that decreasing S has a pronounced effect on 
interface morphology. Dendritic arrays seen in experiments have a certain structure/periodicity to them, that arises 
from underlying crystalline symmetries. For example, in our simulations we observe that the cells constitute a 
hexagonal array. When S is large, away from the boundaries the diffusion field surrounding each cell tip obeys this 
symmetry, and the optimal Ai is selected. As S decreases however, the diffusion field becomes increasingly asymmetric 
due to interaction with the boundaries at x — and x — 6. In particular, solute rejection decreases in the vertical 
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FIG. 9: Comparison of binary alloy simulations with the phase-field model, and experimental data of Liu et al. [2^] . Solid 
symbols correspond to the phase-field model and the open symbols are experimental data. The line denotes a relationship of 
the form pupVup = C between tip radius and velocity, where C can be expressed in terms of process parameters. 



plane x-z, while increasing in the horizontal plane x-y. Increased solute accumulation between cells in x-y, contributes 
to an increase in Ai. However, since Vup is fixed by the pulling speed, and a certain rate of solute rejection needs to 
be maintained, the tips tend to grow sharper as Ai increases. It is precisely this effect that causes the operating state 
to deviate from theoretical predictions. 

To check for what value of S our results deviate from theory, let us assume the steadily growing array at S = 64, 



Vp — 1.5 (Fig. 8(a) I to be one in which cells close to the plane z = 32 are "free" from boundary effects. We estimate 
the c ell sp acing in that plane to be A{ ~ 18, where the superscript 'f denotes "free". When S = 16 and = 1.5 
(Fig. 8(b) I, we notice that our data points start departing from the theoretical prediction, viz. we see a relationship 
of the type PupVup — C, where a > 2. This suggests that agreement with theory deteriorates as S ^ A(, which is 
what Liu et al. concluded from their experiments. A more precise form of the above criterion can be obtained my 
making a careful study of Ai as a function of S, for different Vp, and obtaining a criterion based on a least squares 
fit to the deviation points (as we did with the pure material). Increasing pulling velocity suppresses tip-boundary 
interaction by reducing the thickness of the diffusion boundary layer D/Vup, similar to the effect of melt convection 
in pure materials, and should induce a leftward shift in the curves. 



V. CONCLUDING REMARKS 



We have investigated the role of confinement on solidification in both pure materials and binary alloys. Our 
simulations show that, for equi-axed growth in a pure material, the dendrite's operating state is affected when the 
container dimension S, approaches the scale of the diffusion field D/Vup near the tip. For directionally solidified 
binary alloys, confinement effects become important when 5 is of the order of the primary dendrite spacing Ai. Where 
applicable, one needs to consider the influence of these interactions when comparing experimental data with theoretical 
models that do not account for confinement effects. 

It is notable that we were able to make meaningful comparisons with real experimental results using the phase-field 
model for the alloy. In particular, the agreement obtained in the trends shown by the dendrite tips for different 5 and 
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Vp is very encouraging, and is a testament to the power of phase-field modeling. Some ambiguity remains in classifying 
the computed microstructure as cells or dendrites. Since our results correlated well with dendrite data, we would like 
to think of them as dendrites. Perhaps, incorporating random fluctuations in the phase-field model will resolve this 
issue. We also found tip-splitting inducing oscillations above certain values of the pulling speed. We are unsure as to 
whether this instability has a physical meaning. In experiments, it is seen that increasing Vp causes a decrease in Ai 
and ptip for a steadily growing interface, and our simulations capture this effect, viz. as Vp increases the cells split in 
a manner that produces a more stable configuration with a smaller Ai and pup. However, when Vp > 2.5, it appears 
that an optimal configuration is not possible in our system, given the constraints on its size and boundary conditions. 
We speculate that larger domains should allow for more stable cell configurations at higher pulling speeds. This issue 
too needs further investigation. 

We observed a change in dimensionality of the liquid-solid interface for certain values of S and Vup, when zero 
flux boundary conditions were imposed on the phase-field variable. There is some experimental evidence of this 
phenomenon in the literature. Liu and Kirkaldy [3^ reported a 2-D to 3-D transition in their experiments on a 
SCN-Salol mixture. In their directional solidification experiments in a cell of fixed height {6 = 28pm), they found 
this transition to occur at a driving velocity of 10.8 pm/sec. At a lower driving velocity of 7.6 /im/sec, the dendrites 
looked two dimensional. In our analysis of directional solidification, we found at 5 = 4, the cells underwent a 2-D to 
3-D transition as the pulling velocity was changed from 1 to 1.5. The significance of this result is that through an 
appropriate selection of 6 and Vp in experiments, it should be possible to obtain almost two dimensional dendritic 
arrays in materials that favor wetting. Such experiments will permit more favorable comparisons with 2-D dendrite 
growth theories, since finite boundary effects along the z axis cease to impact the growth. 
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